## pseudo code for ADMB M2
rm(list=ls())
setwd("/Users/kkari/Documents/science/SVN/MSM/svn (trunk)/")
load("pseudoCode.RData")
source("~/Documents/science/SVN/MSM/svn (trunk)/KIRM2/MSMration.R")
outfile<-"/M2.dat"
LWdata2<-read.csv("/Users/kkari/Documents/science/Projects/MSM/msm_data/LengthAtAge2.csv")
LWdata2$LengthCM<-LWdata2$LLENGTH/10
test<-strsplit(as.character(LWdata2$STARTTIME),split="-")
SSTfut<-read.csv("/Users/kkari/Documents/science/Projects/MSM/msm_data/T_future.csv")
Tdev<-SSTfut[,2]-SSTfut[1,2]
rand_Tfut<-matrix(0,100,length(SSTfut[,1]))
rand_Tfut[1,]<-rnorm(length(SSTfut[,1]),mean=mean.na(prey$GearTemp),sd=sd.na(prey$GearTemp)/2)
plot(c(as.numeric(names(meanBT)),SSTfut[,1]),c(meanBT,rand_Tfut[1,]+Tdev),type="l",ylim=c(0,10))
for(i in 2:100){
	rand_Tfut[i,]<-rnorm(length(SSTfut[,1]),mean=mean.na(prey$GearTemp),sd=sd.na(prey$GearTemp)/2)
	lines(c(as.numeric(names(meanBT)),SSTfut[,1]),c(meanBT,rand_Tfut[i,]+Tdev),type="l",col="gray")
}
	lines(c(as.numeric(names(meanBT)),SSTfut[,1]),c(meanBT,apply(rand_Tfut,2,mean.na)+Tdev),type="l",col="black",lwd=2)
	lines(c(as.numeric(names(meanBT)),SSTfut[,1]),c(meanBT,apply(rand_Tfut,2,mean.na)+1.95*apply(rand_Tfut,2,sd.na)+Tdev),type="l",col="black",lwd=2,lty=2)
	lines(c(as.numeric(names(meanBT)),SSTfut[,1]),c(meanBT,apply(rand_Tfut,2,mean.na)-1.95*apply(rand_Tfut,2,sd.na)+Tdev),type="l",col="black",lwd=2,lty=2)
TempC_fut<-apply(rand_Tfut,2,mean.na)+Tdev
Tfut_upper<-apply(rand_Tfut,2,mean.na)+1.95*apply(rand_Tfut,2,sd.na)+Tdev
Tfut_lower<-apply(rand_Tfut,2,mean.na)-1.95*apply(rand_Tfut,2,sd.na)+Tdev

MO<-0
YEAR<-0
MO.DAY<-0
for (i in 1:length(test)){
	MO[i]<-test[[i]][2]
	YEAR[i]<-as.numeric(test[[i]][3])
	MO.DAY[i]<-as.numeric(test[[i]][1])
}
LWdata2$MO<-MO
LWdata2$MO.DAY<-MO.DAY
YEARold<-YEAR
unique(YEAR)
YEAR[YEAR>=82]<-YEAR[YEAR>=82]+1900
YEAR[YEAR<82]<-YEAR[YEAR<82]+2000
years<-sort(unique(YEAR))
LWdata2$YEAR<-YEAR
spp<-unique(LWdata2$GOAPOLL_PRED)

sizes<-seq(0,120,1)
dat<-LWdata2[LWdata2$GOAPOLL_PRED==spp[3],]
plk.N<-hist(dat$LengthCM,breaks=sizes)$counts
dat<-LWdata2[LWdata2$GOAPOLL_PRED==spp[2],]
pcod.N<-hist(dat$LengthCM,breaks=sizes)$counts
dat<-LWdata2[LWdata2$GOAPOLL_PRED==spp[1],]
atf.N<-hist(dat$LengthCM,breaks=sizes)$counts
NatAge<-(data.frame(plk.N,pcod.N,atf.N))
plot(NatAge)
rownames(NatAge)<-sizes
mn.WAplk.jim<-c(0.367,0.520,0.650,0.761,0.871,0.988,1.119,1.207,1.308,1.395,1.448,1.473,1.544)
names(mn.WAplk.jim)<-c(3:14,"15+")
LAplkAsmt<-exp((log(mn.WAplk.jim)-( log(sp.dat$LW.a[1]))/sp.dat$LW.b[1] ) )

#calculate L at Age with ordered logistic regression
spp<-c("WALLEYE POLLOCK","PACIFIC COD","ARROWTOOTH FLOUNDR")
########### Pollock
dat1<-LWdata2[LWdata2$GOAPOLL_PRED==spp[1],]
dat<-na.omit(data.frame(AGE=dat1$AGE,LengthCM=dat1$LengthCM))
ages<-sort(unique(dat$AGE))
datplk<-dat
agesplk<-ages
m.p.plk<-polr(factor(age)~length,data=list(age=dat$AGE,length=dat$LengthCM),method="logistic",Hess = FALSE)  # ordered probit model - proportional odds logistic regression
t.m.p.plk<-data.frame(summary(m.p.plk)[1])
LatA.plk<-cbind(t.m.p.plk[-1,1]/t.m.p.plk[1,1],agesplk[-length(agesplk)])

dat1<-LWdata2[LWdata2$GOAPOLL_PRED==spp[1],]
dat<-na.omit(data.frame(AGE=dat1$AGE,weight=dat1$WEIGHT))
ages<-sort(unique(dat$AGE))
natage.plk<-tapply(dat$AGE,dat$AGE,length.na)
m.p.plk.w<-polr(factor(age)~weight,data=list(age=dat$AGE,weight=dat$weight),method="logistic",Hess = FALSE)  # ordered probit model - proportional odds logistic regression
t.m.p.plk.w<-data.frame(summary(m.p.plk.w)[1])
catt<-strsplit(rownames(t.m.p.plk.w),split="|")
WatA.plk.w<-data.frame(weight=t.m.p.plk.w[-1,1]/t.m.p.plk.w[1,1],age=ages[-1])


########### Pcod
dat1<-LWdata2[LWdata2$GOAPOLL_PRED==spp[2],]
dat<-na.omit(data.frame(AGE=dat1$AGE,LengthCM=dat1$LengthCM))
ages<-sort(unique(dat$AGE))
datpcod<-dat
agespcod<-ages
natage.pcod<-tapply(dat$AGE,dat$AGE,length.na)
m.p.pcod<-polr(factor(age)~length,data=list(age=dat$AGE,length=dat$LengthCM),method="logistic")  # ordered probit model - proportional odds logistic regression
t.m.p.pcod<-data.frame(summary(m.p.pcod)[1],Hess=FALSE)
LatA.pcod<-cbind(t.m.p.pcod[-1,1]/t.m.p.pcod[1,1],agespcod[-length(agespcod)])

########### Arrowtooth
dat1<-LWdata2[LWdata2$GOAPOLL_PRED==spp[3],]
dat<-na.omit(data.frame(AGE=dat1$AGE,LengthCM=dat1$LengthCM))
ages<-sort(unique(dat$AGE))
datatf<-dat
agesatf<-ages
natage.atf<-tapply(dat$AGE,dat$AGE,length.na)
m.p.atf<-polr(factor(age)~length,data=list(age=dat$AGE,length=dat$LengthCM),method="logistic",Hess = FALSE)  # ordered probit model - proportional odds logistic regression
t.m.p.atf<-data.frame(summary(m.p.atf)[1])
LatA.atf<-cbind(t.m.p.atf[-1,1]/t.m.p.atf[1,1],agesatf[-length(agesatf)])

## expand to reach 21 ages
ages<-1:21
LatA.atf2<-data.frame(L=0,Age=ages)
LatA.atf2[na.omit(match(ages,LatA.atf[,2])),1]<-LatA.atf[na.omit(match(LatA.atf[,2],ages)),1]
LatA.atf2[18:21,1]<-c(82,84,86,88)
pcodAsmt<-data.frame(age=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20),
Length=c(23.08, 37.44, 48.9, 58.04, 65.33, 71.15, 75.79, 79.49, 82.44, 84.79, 86.67, 88.16, 89.36, 90.31, 91.07, 91.68, 92.16, 92.54, 92.85, 93.3))
## plot results

graphics.off()
par(mfrow=c(3,1))
plot(datplk$AGE,datplk$LengthCM, pch=16,cex=.8, main="pollock")
points(LatA.plk[,2],LatA.plk[,1],col="red",pch=16)
lines(3:15,LAplkAsmt*10,col="red")
plot(datpcod$AGE,datpcod$LengthCM, pch=16,cex=.8, main="pcod")
points(LatA.pcod[,2],LatA.pcod[,1],col="red",pch=16)
lines(pcodAsmt,col="red")
plot(datatf$AGE,datatf$LengthCM, pch=16,cex=.8, main="arrowtooth")
points(LatA.atf[,2],LatA.atf[,1],col="red",pch=16)



i<-1
	dat<-LWdata2[LWdata2$GOAPOLL_PRED==spp[i],]
	plot(dat$AGE,dat$WEIGHT)
	points(LatA.plk[,2],sp.dat$LW.a[i]*(LatA.plk[,1]^sp.dat$LW.b[i]),col="red",pch=16)
	cbind(LatA.plk,weightKG=(sp.dat$LW.a[i]*(LatA.plk[,1]^sp.dat$LW.b[i]))/1000)

i<-2
	dat<-LWdata2[LWdata2$GOAPOLL_PRED==spp[i],]
	plot(dat$AGE,dat$WEIGHT)
	points(LatA.pcod[,2],sp.dat$LW.a[i]*(LatA.pcod[,1]^sp.dat$LW.b[i]),col="red",pch=16)
i<-3
	dat<-LWdata2[LWdata2$GOAPOLL_PRED==spp[i],]
	plot(dat$AGE,dat$WEIGHT)
	points(LatA.atf[,2],sp.dat$LW.a[i]*(LatA.atf[,1]^sp.dat$LW.b[i]),col="red",pch=16)
#########################################################
## 2. Calculate Upaij observed for each pred - eg. test.21$Upaij is the Upaij for pcod eating pollock
#########################################################
prey.l<-seq(1,120,1)
pred.l<-seq(1,120,1)
nprey<-2
npred<-3
spp.short<-c("plk","pcod","atf")
Kobs<-array(0,c(npred,nprey,length(pred.l)))
Utmp<-matrix(0,length(pred.l),length(prey.l))
colnames(Utmp)<-paste("prey ",prey.l,sep="")
rownames(Utmp)<-paste("pred ",pred.l,sep="")
preysppl<-c("W. Pollock","P. Cod","Arrowtooth")
preyspp<-c("WALLEYE.POLLOCK_Wt","PACIFIC.COD_Wt","ARROWTOOTH_Wt")
Uobs<-array(NA,c(npred,nprey,length(pred.l),length(prey.l)))# U(predd,preyy
par(mfrow=c(3,1))

for (predd in 1:npred){	
	eval(parse(text=paste("sub.datl<-",sp.dat$spp[predd],".l",sep="")))
	sub.datl$Wt<-sp.dat$LW.a[predd]*(sub.datl$PRED_LEN^sp.dat$LW.b[predd])
	eval(parse(text=paste("sub.dat<-",sp.dat$spp[predd],sep="")))
	sub.dat$Wt<-sp.dat$LW.a[predd]*(sub.dat$Length^sp.dat$LW.b[predd])
	for(preyy in 1:nprey){
		Utmp<-matrix(0,length(pred.l),length(prey.l))
		prey.lwt<-sp.dat$LW.a[preyy]*(prey.l^sp.dat$LW.b[preyy])
		colnames(Utmp)<-paste(spp.short[preyy],prey.l,sep=" ")
		rownames(Utmp)<-paste(spp.short[predd],pred.l,sep=" ")
		sub.datl.p<-sub.datl[sub.datl$ECOPATH_PREY==preysppl[preyy],]
		sub.dat.p<-sub.dat[sub.dat$ECOPATH_PREY==preyspp[preyy],]
		Ktmp<-eval(parse(text=paste("tapply(sub.dat$",preyspp[preyy],"/sub.dat$TotWt,sub.dat$Length,mean.na)",sep="")))
		Kobs[predd,preyy,na.omit(match(as.numeric(names(Ktmp)),pred.l))]<-Ktmp[is.na(match(as.numeric(names(Ktmp)),pred.l))==FALSE]
		for (l in 1:length(pred.l)){
			subb<-sub.datl.p[sub.datl.p$PRED_LEN==pred.l[l],]
			tt<-tapply(subb$FREQ,subb$PREY_SIZE_CM,sum.na) # frequency in diet
			tt<-tt*sp.dat$LW.a[preyy]*(as.numeric(names(tt))^sp.dat$LW.b[preyy]) # frequency by mass
			tt<-tt/sum(tt) # proportion of total in diet by weight
			Utmp[l,na.omit(match(as.numeric(names(tt)),prey.l))]<-tt[is.na(match(as.numeric(names(tt)),prey.l))==FALSE]
			Utmp[l,]<-Utmp[l,]*Kobs[predd,preyy,l]  # convert to proportion by weight
			rm(tt)
			rm(subb)
		}
		Uobs[predd,preyy,,]<-Utmp
	}
}



## plot of size preference curve
X11(width=5,height=3)
#layout(rbind(
#		c(1,0,4),
#		c(2,0,5),
#		c(3,0,6))
#		,widths=c(1,0.01,1),heights=rep(1,3))
lim<-60
l1<-4
l2<-30
for (i in 1:length(x)){if(x[i]>lim){y[i]<-exp(-((x[i]-lim)^2)/l1)}else{y[i]<-exp(-((lim-x[i])/l2)^2)}}
plot(x,y,type="l",lwd=2, ylab="preference",xlab="prey length",axes=FALSE)
abline(v=lim,lwd=2,lty=2)
axis(1);axis(2)
lim<-30
for (i in 1:length(x)){if(x[i]>lim){y[i]<-exp(-((x[i]-lim)^2)/l1)}else{y[i]<-exp(-((lim-x[i])/l2)^2)}}
lines(x,y,type="l",lwd=2, col="red")
abline(v=lim,lwd=2,lty=2,col="red")
lim<-120
l1<-10
l2<-50
for (i in 1:length(x)){if(x[i]>lim){y[i]<-exp(-((x[i]-lim)^2)/l1)}else{y[i]<-exp(-((lim-x[i])/l2)^2)}}
lines(x,y,type="l",lwd=2, col="blue")
abline(v=lim,lwd=2,lty=2,col="blue")


source("~/Documents/science/R_funKir/contplot.R")
#__
#library(mgcv)
colss<-c(rep("white",2),heat.colors(20)[20:0],rep("red",20))
graphics.off()
Uhat<-array(0,c(npred,nprey,length(pred.l),length(prey.l)))# U(predd,preyy
Uhatcum<-array(0,c(nprey,length(pred.l),length(prey.l)))# U(predd,preyy
Udatcum<-Uhatcum
klim<-c(80,110,78)
#klim<-c(150,150,150)

for (predd in 1:npred){
	for (preyy in 1:1){
		eval(parse(text=paste("sub.datl<-",sp.dat$spp[preyy],".l",sep="")))
		eval(parse(text=paste("sub.dat<-",sp.dat$spp[preyy],sep="")))

		redd1<-.4
		redd2<-.5
		Udat<-Uobs[predd,preyy,,]
		subb<-tapply(sub.datl$PRED_LEN,sub.datl$PREDJOIN,max.na)
		ll<-tapply(subb,subb,length.na)
		ll.w<-tapply(sub.dat$Length,sub.dat$Length,length.na)
		#ll<-ll/sum(ll)
		prey.lwt<-sp.dat$LW.a[preyy]*(prey.l^sp.dat$LW.b[preyy])
		ll2<-colSums(Udat)/sum(colSums(Udat))
		envL<-rep(0,length(ll2))
		envL2<-envL
		envL3<-envL	
		envL[match(names(ll),prey.l)]<-ll
		envLW<-envL*prey.lwt
		envLW<-envLW/sum(envLW)
		envL<-envL/sum(envL)

		envL2[match(as.numeric(names(ll.w)),prey.l)]<-ll.w
		envLW2<-envL2*prey.lwt
		envLW2<-envLW2/sum(envLW2)
		envL2<-envL2/sum(envL2)
		subLW<-LWdata2[LWdata2$GOAPOLL_PRED==spp[preyy],]
		subLW$LengthCM<-round(subLW$LengthCM)
		ll3<-tapply(subLW$LengthCM,subLW$LengthCM,length.na)
		envL3[na.omit(match(as.numeric(names(ll3)),prey.l))]<-ll3
		envLW3<-envL3*prey.lwt
		envLW3<-envLW3/sum(envLW3)
		envL3<-envL3/sum(envL3)
		
		pref<-matrix(0,length(pred.l),length(prey.l))
		switcha<-pref
		lim<-sp.dat$lim.glm.a[predd]+sp.dat$lim.glm.b[predd]*pred.l
		lim[lim<0]<-0
		Uhat.tmp<-pref
		l1<-4
		l2<-200
		l3<-c(.55,.55,1)# size pref
		Kobs.g<-Kobs
		Kobs.g[Kobs==0]<-0.00001
		Kobs.g[Kobs==1]<-0.99999
		rr<-which(colnames(sub.dat)==sp.dat$spp.p.prey[preyy])
		sub.sub.dat<-sub.dat[sub.dat$Length<=klim[predd],]
		K.sub<-sub.sub.dat[,rr]/sub.sub.dat$TotWt
		K.sub[K.sub==1]<-0.99999
		K.sub[K.sub==0]<-0.00001
		gam.dat<-na.omit(data.frame(K=K.sub,L=sub.sub.dat$Length))
#		Kglm<-gam(K~s(L),data=data.frame(K=Kobs.g[predd,preyy,][pred.l<=klim[predd]],L=pred.l[pred.l<=klim[predd]]),family=Gamma(link="logit"))
#		Kgam<-gam(K~s(L),data=gam.dat,family=Gamma(link="logit"))

#		Khat<-predict.gam(Kgam,type="response", newdata=data.frame(L=pred.l))
		Khat<-exp(Kglm_a[predd,preyy]+((pred.l)*Kglm_b[predd,preyy]))
		Khat[Khat<0]<-0
		Khat[Khat>1]<-1
		Khat[pred.l>klim[predd]]<-0
		for (j in 1:length(pred.l)){
			part1<-rep(0,length(prey.l))
			for (a in 1:length(prey.l)){
				if(a>lim[j]){
					switcha[j,a]<-exp(-((max(0,a-lim[j]))^1)/l1)
					
					#
				}else{
					switcha[j,a]<-1# exp(-((max(0,lim[j]-a))^2)/l2)#1-0.5*((lim[j]-a)/lim[j]) # 
					
				}
				pref[j,a]<-exp(-((a-(lim[j]*l3[predd]))^2)/l2) #exp(-(a-(j*l3[predd])^2)/l2) #
				part1[a]<-pref[j,a]*envL3[a]*switcha[j,a]
			}
			
			Uhat.tmp[j,]<-Khat[j]*part1/sum(part1)
		}
		
		Uhat[predd,preyy,,]<-Uhat.tmp
		Uhatcum[preyy,,]<-Uhatcum[preyy,,]+Uhat.tmp
		Udatcum[preyy,,]<-Udatcum[preyy,,]+Udat
		
		cont.plot(main="Upaij",predname=preysppl[predd],preyname=preysppl[preyy],nlev=10,
			redd=redd1,zdat=Udat,zdat.hat=Uhat.tmp,xdat=pred.l,
			ydat=prey.l,yprey=cbind(apply(Udat,2,sum.na),envL3,envLW3,apply(Uhat.tmp,2,sum.na)),
			xprey=prey.l,xpred=pred.l,ypred=cbind(apply(Udat,1,sum.na),apply(Uhat.tmp,1,sum.na)),
			xlimm=c(0,120),ylimm=c(0,120),zlimm=c(0,.1),splts=100,adj=4.5,padj=2.8,lim=lim)

	cont.plot(main="Upaij_cum",predname=preysppl[predd],preyname=preysppl[preyy],nlev=10,
			redd=redd1,zdat=Udatcum[preyy,,],zdat.hat=Uhatcum[preyy,,],xdat=pred.l,
			ydat=prey.l,yprey=cbind(apply(Udatcum[preyy,,],2,sum.na),envL3,envLW3,apply(Uhatcum[preyy,,],2,sum.na)),
			xprey=prey.l,xpred=pred.l,ypred=cbind(apply(Udatcum[preyy,,],1,sum.na),apply(Uhatcum[preyy,,],1,sum.na)),
			xlimm=c(0,120),ylimm=c(0,120),zlimm=c(0,.1),splts=100,adj=4.5,padj=2.8,lim=lim)

	}
	
}
a<-1:100
lim<-50
tt<-rep(0,length(a))
for (i in 1:length(a)){

tt[i]<-exp(-((max(0,a[i]-lim))^2)/4)
}

Uobs[is.na(Uobs)]<-0
Kobs[is.na(Kobs)]<-0
# create srv_age_sizes for pollock
# create srv_age_err for pollock

means.plk<-tapply(datplk$LengthCM,datplk$AGE,mean)
sd.plk<-tapply(datplk$LengthCM,datplk$AGE,sd)
bins<-round(means.plk)[-1]
bins<-bins[1:12]
#bins<-LatA.plk[-1,1]
tt<-matrix(0,12,length(bins))
aa<-hist(datplk$LengthCM[datplk$AGE==1],plot=FALSE)
plot(aa$breaks[-length(aa$density)],aa$density,type="l",ylim=c(0,.2),xlim=c(0,100))
plot(density(datplk$LengthCM[datplk$AGE==1]),ylim=c(0,.2),xlim=c(0,100))
lines(density(rnorm(1000,mean=means.plk[1],sd=sd.plk[1])),col="red")
aa<-density(datplk$LengthCM[datplk$AGE==1])
for (a in 1:12){
	aa<-hist(datplk$LengthCM[datplk$AGE==a],plot=FALSE)
	lines(density(datplk$LengthCM[datplk$AGE==a]),ylim=c(0,.2),xlim=c(0,100))

	##lines(aa$breaks[-length(aa$density)],aa$density,type="l")
	tt[a,]<-dnorm(bins,mean=means.plk[a+1],sd=sd.plk[a+1])
	tt[a,]<-tt[a,]/sum(tt[a,])
	lines(density(rnorm(1000,mean=means.plk[a+1],sd=sd.plk[a+1])),col="red")
}


############################################################
### output to .dat file 
############################################################
outfile<-"/Users/kkari/Documents/science/SVN/MSM_copy/stomach.dat"
onn<-1
	nages<-c(12,12,21)
	nyrs<-29
if(onn==1){
	app1<-FALSE
	cat("#Kirs : run bioenergetics or not ",file=outfile,append=app1,sep="\n")
	cat(1,file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#npred2 : number of predator spp ",file=outfile,append=app1,sep="\n")
	cat(3,file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#useWt : if ==1 assign relative prop of prey in the diet according to relative biomass in the system.,otherwise the model will use size freq",file=outfile,append=TRUE,sep="\n")
	cat(c(1,1,1),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#C_model :if ==1, the use Cmax*fT*P, if !=1 then use livingston ration approach ",file=outfile,append=TRUE,sep="\n")
	cat(c(1,1,1),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#nprey2: number of prey spp ",file=outfile,append=TRUE,sep="\n")
	cat(2,file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#nspp2 : number of species ",file=outfile,append=TRUE,sep="\n")
	cat(3,file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#nyrs2 : number of years ",file=outfile,append=TRUE,sep="\n")
	cat(32,file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#nages2 : number of unique age classes for each species ",file=outfile,append=TRUE,sep="\n")
	cat(nages,file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
##CMAX paramters
	cat("#Pvalue : only used if C_model ==1 , proportion of Cmax; Pvalue is P in Cmax*fT*P",file=outfile,append=TRUE,sep="\n")
	cat(c(1,.3,1),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#Ceq : which Comsumption equation to use",file=outfile,append=TRUE,sep="\n")
	cat(model.parms$Ceq[c(1,4,5)],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#CA : Wt specific intercept of Cmax=CA*W^CB",file=outfile,append=TRUE,sep="\n")
	cat(model.parms$CA[c(1,4,5)],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#CB : Wt specific slope of Cmax=CA*W^CB",file=outfile,append=TRUE,sep="\n")
	cat(model.parms$CB[c(1,4,5)],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#Qc : used in fT, QC",file=outfile,append=TRUE,sep="\n")
	cat(model.parms$CQ[c(1,4,5)],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#Tco : used in fT, thermal optimum",file=outfile,append=TRUE,sep="\n")
	cat(model.parms$Tco[c(1,4,5)],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#Tcm : used in fT thermal max",file=outfile,append=TRUE,sep="\n")
	cat(model.parms$Tcm[c(1,4,5)],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#Tcl : used in fT eq 3 thermallimit",file=outfile,append=TRUE,sep="\n")
	cat(c(30,30,30),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#CK1 : used in fT eq 3, temp where C is .98 max (ascending)",file=outfile,append=TRUE,sep="\n")
	cat(c(5,5,5),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")	
	cat("#CK4 : used in fT eq 3, limit where C is .98 max (decending)",file=outfile,append=TRUE,sep="\n")
	cat(c(17,17,17),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")

##END CMAX paramters
#	cat("#nlengths : number of unique age classes for each species ",file=outfile,append=TRUE,sep="\n")
#	cat(length(sizes[-1]),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")	
#	cat("#lengths : length classess for AvgN length and number datafile - based on average values for all years ",file=outfile,append=TRUE,sep="\n")
#	cat(sizes[-1],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
#	cat("#avgN : mean number in each classess for AvgN length and number datafile - based on average values for all years ",file=outfile,append=TRUE,sep="\n")
#	cat("#pollock ",file=outfile,append=TRUE,sep="\n")
#	for(i in 1:32){cat(NatAge[,1],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")}
#	cat("#pcod ",file=outfile,append=TRUE,sep="\n")
#	for(i in 1:32){cat(NatAge[,2],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")}
#	cat("#arrowtooth  ",file=outfile,append=TRUE,sep="\n")
#	for(i in 1:32){cat(NatAge[,3],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")}
	cat("#Kglm_a : LK a regression coef for K=exp(a[pred,prey]+L*b[pred,prey]), returns a vector with 1xnsize, rows=pred, cols=prey ",file=outfile,append=TRUE,sep="\n")
	for(i in 1:3){cat((sp.dat$a.glm[i,]),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")}
	cat("#Kglm_b : LK b regression coef for K=exp(a[pred,prey]+L*b[pred,prey]), rows=pred, cols=prey ",file=outfile,append=TRUE,sep="\n")
	for(i in 1:3){cat((sp.dat$b.glm[i,]),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")}
	#temp
	cat("#nTyrs : number Temperature years ",file=outfile,append=TRUE,sep="\n")
	cat(length(t(Temp.data$Year)),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#Tyrs : Temperature years ",file=outfile,append=TRUE,sep="\n")
	cat(t(Temp.data$Year),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#TempC : Temperature ",file=outfile,append=TRUE,sep="\n")
	cat(t(Temp.data$EBS_bottomT),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#S_a : a,L,L^2,L^3,L^4,L^5 (rows)coef for mean S=a+b*L+b2*L*L, whith a cap at 80cm for each pred spp(cols)",file=outfile,append=TRUE,sep="\n")
	cat((sp.dat$aL),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat((sp.dat$bL),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat((sp.dat$bL2),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat((sp.dat$bL3),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat((sp.dat$bL4),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat((sp.dat$bL5),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	#limit
	cat("#aLim and bLim : a&b regression coef for W=a*L^b ",file=outfile,append=TRUE,sep="\n")
	cat((sp.dat$lim.glm.a),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat((sp.dat$lim.glm.b),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	#diet preference
	#LW regression coef for each spp
	cat("#aLW : LW a&b regression coef for W=a*L^b ",file=outfile,append=TRUE,sep="\n")
	cat((sp.dat$LW.a),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat((sp.dat$LW.b),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#maxK : nsppxnspp matrix of maximum diet proportions for each predd,preyy combo (used for broken stick expon)",file=outfile,append=TRUE,sep="\n")
	for(i in 1:3){cat(round(sp.dat$maxK[i,],2),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")}
	cat("#Length at Age matrix for each species  ",file=outfile,append=TRUE,sep="\n")
	cat("#LA_age :",file=outfile,append=TRUE,sep="\n")
	cat(1:nages[1],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat(1:nages[2],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat(1:nages[3],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#LA_Lengths : rounded to the nearest integer",file=outfile,append=TRUE,sep="\n")
	cat(round(LatA.plk[,1][match(1:nages[1],LatA.plk[,2])]),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat(round(LatA.pcod[,1][match(1:nages[2],LatA.pcod[,2])]),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat(round(LatA.atf2[,1][match(1:nages[3],LatA.atf2[,2])]),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#nyrs_fut : to make sure file reads data correctly",file=outfile,append=TRUE,sep="\n")
	cat(length(TempC_fut),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#TempC_fut_yrs : to make sure file reads data correctly",file=outfile,append=TRUE,sep="\n")
	cat(SSTfut[,1],file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#TempC_fut : mean, upper and lower to make sure file reads data correctly",file=outfile,append=TRUE,sep="\n")
	cat(TempC_fut,file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat(Tfut_upper,file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat(Tfut_lower,file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
	cat("#test_stomach : to make sure file reads data correctly",file=outfile,append=TRUE,sep="\n")
	cat(12345,file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")


	## output observed U
	cat("#Uobs",file=outfile,append=TRUE,sep="\n")
	for(predd in 1:npred){
		cat(paste("#Uobs","_",spp.short[predd],sep=""),file=outfile,append=TRUE,sep="\n")
		for (preyy in 1:nprey){
			Udat<-Uobs[predd,preyy,,]
			for (i in 1: length(Udat[,1])){
				cat(round(Udat[i,],4),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
			}
		}
	}
	## output observed K
	cat("#Kobs",file=outfile,append=TRUE,sep="\n")
	for(predd in 1:npred){
		Kdat<-Kobs[predd,,]
		for (i in 1: length(Kdat[,1])){
			cat(round(Udat[i,],4),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")
		}
	}
	cat("#srv_age_size :plk to make sure file reads data correctly",file=outfile,append=TRUE,sep="\n")
	cat(bins,file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")

	cat("#srv_age_err :plk to make sure file reads data correctly",file=outfile,append=TRUE,sep="\n")
	for (i in 1:12){cat(round(tt[i,],4),file=outfile,append=TRUE,sep="\t");cat("",file=outfile,append=TRUE,sep="\n")}

}
############################################################
BELOW NEEDS TO BE UPDATED!
############################################################


###Pseudo code for M2
Kglm_a<-sp.dat$a.glm
Kglm_b<-sp.dat$b.glm
Tyrs<-Temp.data$Year;TempC<-Temp.data$EBS_bottomT
# mean stom wt (as fraction of body weight)
stomLength<-stomLength
#StomWt<-StomWt  # for now read it in as data
S_a<-sp.dat$aL
S_b<-sp.dat$bL
S_b2<-sp.dat$bL2

aLim<-sp.dat$lim.glm.a
bLim<-sp.dat$lim.glm.b

nprey<-2 ## though zeros across the board for N==3
npred<-3
pred_size<-1:140
prey_size<-1:160
n_pred_size<-length(pred_size)
n_prey_size<-length(prey_size)

# size frequency
# limit
# diet preference
# LW regression coef for each spp
aLW<-sp.dat$LW.a
bLW<-sp.dat$LW.b
maxK<-maxK

##### !!!!!! NEED TO ADD SIZE CONVERSION MATRIX
age_size<-rbind(tapply(pollock$Length,pollock$pred.age,min.na)[-1],tapply(pcod$Length,pcod$pred.age,min.na)[-1])

age_size[1,which(is.na(age_size[1,]))]<-age_size[1,which(is.na(age_size[1,]))-1]+2 # replace missing values
age_size[2,which(is.na(age_size[2,]))]<-age_size[2,which(is.na(age_size[2,]))-1]+2 # replace missing values
age_size<-rbind(age_size,age_size[2,])
######

skipp<-0
if (skipp==0){
	## W at Age for each spp:
	##POLLOCK
	dat1<-data.frame(AGE=plk_WA$AGE,WEIGHT=plk_WA$WEIGHT,LENGTH=plk_WA$LENGTH)
	dat1<-na.omit(dat1)
	ages<-sort(unique(dat1$AGE))
	m.p<-polr(factor(age)~weight,data=list(age=dat1$AGE,weight=dat1$WEIGHT),method="logistic")  # ordered probit model - proportional odds logistic regression
	t.m.p<-data.frame(summary(m.p)[1])
	graphics.off()
	par(mfrow=c(3,1))
	plot(dat1$AGE,dat1$WEIGHT,pch=16,cex=.8, main="pollock")
	points(ages[-1],t.m.p[-1,1]/t.m.p[1,1],type="b",pch=16,col="red",lwd=1,main="pollock")
	WatA.plk<-data.frame(Weight=t.m.p[-1,1]/t.m.p[1,1],Age=ages[-1])
	
	dat1<-data.frame(AGE=plk_WA$AGE,WEIGHT=plk_WA$WEIGHT,LENGTH=plk_WA$LENGTH)
	ages<-sort(unique(dat1$AGE))
	m.p<-polr(factor(age)~length,data=list(age=dat1$AGE,length=dat1$LENGTH),method="logistic")  # ordered probit model - proportional odds logistic regression
	t.m.p<-data.frame(summary(m.p)[1])
	LatA.plk<-cbind(t.m.p[-1,1]/(10*t.m.p[1,1]),ages[-1])
	WatA.plk$Weights2<-(sp.dat$LW.a[1]*(LatA.plk[,1])^sp.dat$LW.b[1])[-1]
	points(WatA.plk$Age,WatA.plk$Weights2,type="b",pch=16,col="blue",lwd=1,main="pollock")
	##PCOD
	dat1<-data.frame(AGE=pcod_WA$AGE,WEIGHT=pcod_WA$WEIGHT,LENGTH=pcod_WA$LENGTH)
	dat1<-na.omit(dat1)
	ages<-sort(unique(dat1$AGE))
	m.p<-polr(factor(age)~weight,data=list(age=dat1$AGE,weight=dat1$WEIGHT),method="logistic")  # ordered probit model - proportional odds logistic regression
	t.m.p<-data.frame(summary(m.p)[1])
	plot(dat1$AGE,dat1$WEIGHT,pch=16,cex=.8, main="pollock")
	points(ages[-1],t.m.p[-1,1]/t.m.p[1,1],type="b",pch=16,col="red",lwd=1,main="pollock")
	WatA.pcod<-data.frame(Weight=t.m.p[-1,1]/t.m.p[1,1],Age=ages[-1])
	
	##ATF
	dat<-LWdata2[LWdata2$GOAPOLL_PRED=="ARROWTOOTH FLOUNDR",]
	dat1<-data.frame(AGE=dat$AGE,WEIGHT=dat$WEIGHT,LENGTH=dat$LENGTH)
	dat1<-na.omit(dat1)
	ages<-sort(unique(dat1$AGE))
	m.p<-polr(factor(age)~weight,data=list(age=dat1$AGE,weight=dat1$WEIGHT),method="logistic")  # ordered probit model - proportional odds logistic regression
	t.m.p<-data.frame(summary(m.p)[1])
	plot(dat1$AGE,dat1$WEIGHT,pch=16,cex=.8, main="pollock")
	points(ages[-1],t.m.p[-1,1]/t.m.p[1,1],type="b",pch=16,col="red",lwd=1,main="pollock")
	WatA.atf<-data.frame(Weight=t.m.p[-1,1]/t.m.p[1,1],Age=ages[-1])
}

#dummy data
#toy data
nyrs<-10
srv_age_Wt<-WatA.plk[,2]
nages<-c(12,12,21)#rep(dim(age_size)[2],3)
srv_age_Wt<-WatA.plk[,1]
AvgN<-rbind(tapply(pollock$pred.age,pollock$pred.age,length.na)[-1],tapply(pcod$pred.age,pcod$pred.age,length.na)[-1],tapply(pcod$pred.age,pcod$pred.age,length.na)[-1])
AvgN[is.na(AvgN)]<-0

N<-array(0,c(nyrs,nages[1],3))  #pre-allocate dummy N data
biomass<-array(0,c(nyrs,nages[1],3))  #pre-allocate dummy biomass data
for (spp in 1:3){
	for (i in 1:nyrs){
		N[i,,spp]<-rnorm(nages[1],0,.2)+AvgN[spp,]
		biomass[i,,spp]<-N[i,,spp]*aLW[spp]*age_size[spp,]^bLW[spp]
	}
}


srv_age_obs<-matrix(0,npred,nages[1]) # min L A conversion matrix
srv_age_obs[1:2,]<-rbind(tapply(pollock$Length,pollock$pred.age,min.na)[-1],tapply(pcod$Length,pcod$pred.age,min.na)[-1])

srv_age_obs[1,which(is.na(srv_age_obs[1,]))]<-srv_age_obs[1,which(is.na(srv_age_obs[1,]))-1]+2 # replace missing values
srv_age_obs[2,which(is.na(srv_age_obs[2,]))]<-srv_age_obs[2,which(is.na(srv_age_obs[2,]))-1]+2 # replace missing values
srv_age_obs[3,]<-srv_age_obs[2,]





##Pseudo code
#rm(sp.dat)
####################################################
## Local calcs
## Calculate S: mean stomach weight (as prop of body weight) as a function of length 
####################################################
S<-matrix(0,npred,nages[1])  #pre-allocate
Wt<-S
for (predd in 1:npred){
	S[predd,]<-S_a[predd]+S_b[predd]*age_size[predd,]+S_b2[predd]*(age_size[predd,]^2)
		for(j in 1:nages){
			if (S[predd,age_s
	S[predd,age_size[predd,]>80]<-S_a[predd]+S_b[predd]*80+S_b2[predd]*(80^2)
	Wt[predd,]<-aLW[predd]*age_size[predd,]^bLW[predd]
}

####################################################
## Calculate K (diet pref) for each size of predator - can be changed from sizes to ages using median or mean size per age
####################################################
K<-array(0,c(nprey,nages[1],npred)) #pre-allocate


for(predd in 1:npred){
	Ktmp<-matrix(0,nprey,nages[predd])#pre-allocate
	for(preyy in 1:nprey){		
		for(j in 1:nages[predd]){
			Ktmp[preyy,j]<-min(1,exp(Kglm_a[predd,preyy]+(age_size[predd,j]*Kglm_b[predd,preyy])))
		}
		Ktmp[preyy,Ktmp[preyy,]>maxK[predd,preyy]]<-maxK[predd,preyy]
		if(predd==1){Ktmp[preyy,age_size[predd,]>90]<-0}
		Ktmp[,colSums(Ktmp)>1]<-Ktmp[,colSums(Ktmp)>1]/colSums(Ktmp[,colSums(Ktmp)>1])	#re-standardize
	}# end prey	
	K[,,predd]<-Ktmp
}#end pred

####################################################
## Calculate lim upper size limit that each predator can eat (not prey dependent) - can be changed from sizes to ages using median or mean size per age
####################################################
Limit<-matrix(0,npred,nages[1])#pre-allocate

for(predd in 1:npred){
	Limit[predd,]<-aLim[predd]+bLim[predd]*age_size[predd,]
}#end pred

####################################################
##based on AGE Calculate (U) effective prey size freq (Npa/sum(N))*switch -- rows=prey, cols=preds, can be in calcs section if Npa/Np doesn't change each year
####################################################
U<-array(0,c(nages[1],nages[1],nprey,npred)) #pre-allocate
 
for (predd in 1:npred){
	for(preyy in 1:nprey){
		Utmp<-matrix(0,nages[preyy],nages[preyy])
		prey.prop<-AvgN[preyy,]/sum(AvgN[preyy,])##### <- change this to numerical abundance by year !!!!
		for(j in 1:nages[predd]){
			switch1<-rep(1,nages[preyy])
			switch1[age_size[preyy,]>=Limit[predd,j]]<-exp(-(age_size[preyy,][age_size[preyy,]>=Limit[predd,j]]-Limit[predd,j])/4)
			prey.use<-prey.prop*switch1/sum(prey.prop*switch1)
			Utmp[,j]<-prey.use*K[preyy,j,predd]  # can't get to work as matrix mult - need more coffee - prey.use is 140x1, and K is 1x140 =140x140
		}
		U[,,preyy,predd]<-Utmp	
		#eval(parse(text=paste("Upaij.hat.",predd,preyy,"<-Upaij.hat",sep="")))	
	}
}
####################################################
## Calculate annual ration C=24*0.0134*exp(0.0115*TempC)*91.25*S  units are g of food/g of predator
####################################################
C<-array(0,c(nyrs,nages[1],npred)) #pre-allocate
for(i in 1:nyrs){
	for(predd in 1:npred){
		C[i,,predd]=24*0.0134*exp(0.0115*TempC[i])*91.25*S[predd,]
	}
}

####################################################
## Calculate grams of prey eaten (E) units are g of prey E= U*C*N(pred,pred_size)*wt(pred,pred_size)
####################################################
E<-array(0,c(nyrs,nages[1],nprey,npred)) #pre-allocate
Esum<-array(0,c(nyrs,nages[1],nprey)) #pre-allocate
for(i in 2:nyrs){
	for(preyy in 1:nprey){
		Etmp<-matrix(0,npred,nages[1])
		for (predd in 1:npred){
			Etmp[predd,]<-rowSums(t(t(U[,,preyy,predd])*C[i,,predd]*biomass[i-1,,predd]))  ### Note that the array must be flipped to mult and flipped again to get back to prey in rows, pred in cols
		}
		E[i,,preyy,]<-Etmp
	Esum[i,,preyy]<-colSums(Etmp)  # returns a matrix nyrs x n_prey_size
	}
}


####################################################
## LASTLY!! Calculate M2 as E/N
####################################################
M2<-array(0,c(nyrs,nages[1],nprey)) #pre-allocate

for(i in 2:nyrs){
	for(preyy in 1:nprey){
		M2[i,,preyy]<-E[i,,preyy]/(biomass[i-1,,preyy])  # need to make biomass updated in each year....
	}
}

U.hat<-U
####################################################
## Plot U and Uhat
####################################################
# to calculate U obs
# 1. from prey.l data determine the # of prey size j in diet of pred size k
# 2. Convert (1) to Wt of prey size j in diet of pred. size k by multiplying the length (j) by the LW wt relationship W=a*L^b
# 3. Standardize 2 to the colsums of pred size k (such that Wt is Wt/sum(Wt) for each pred size k)
# 4. Multiply (3) by the mean wt of the prey in the diet of pred size k (diet pref), using data from prey (weight data)

#pcod Eat plk
dat<-pcod.l[pcod.l$ECOPATH_PREY=="W. Pollock",]

# 1. from prey.l data determine the # of prey size j in diet of pred size k
	U.tmpold<-pivot.table(dat$PREY_SIZE,dat$PREY_SIZE_CM,dat$PRED_LEN,sum.na,missing=0)
	levelplot(t(U.tmpold))
# 2. Convert (1) to Wt of prey size j in diet of pred. size k by multiplying the length (j) by the LW wt relationship W=a*L^b
	U.tmp<-pivot.table(sp.dat$LW.a[1]*(dat$PREY_SIZE^sp.dat$LW.b[1]),dat$PREY_SIZE_CM,dat$PRED_LEN,sum.na,missing=0)
	levelplot(t(U.tmp))
# 3. Standardize 2 to the colsums of pred size k (such that Wt is Wt/sum(Wt) for each pred size k)
	tt<-colSums(U.tmp)
	for (r in 1:length(U.tmp[,1])){
		U.tmp[r,]<-U.tmp[r,]/tt
	}
	levelplot(t(U.tmp))
# 4. Multiply (3) by the mean wt of the prey in the diet of pred size k (diet pref), using data from prey (weight data)
	kk<-tapply(pcod$WALLEYE.POLLOCK_Wt/pcod$TotWt,pcod$Length,mean.na)
	rr<-colnames(U.tmp)
	cc<-match(rr,names(kk))
	for (r in 1:length(U.tmp[,1])){
		U.tmp[r,]<-U.tmp[r,]*kk[cc]
	}	
	levelplot(t(U.tmp))
	#plot(rowSums(U.tmp), type="l")
U<-matrix(0,100,120)
rownames(U)<-1:100; colnames(U)<-1:120
for (preyL in 1:100){# for each row
	if(any(rownames(U.tmp)==preyL)){
		rr<-which(rownames(U.tmp)==preyL)
		cc<-na.omit(match(colnames(U.tmp),colnames(U)))
		cc2<-na.omit(match(colnames(U),colnames(U.tmp)))
		U[preyL,cc]<-U.tmp[rr,cc2]
		rm(rr)
	}
}
	levelplot(t(U))
	
	levelplot(t(U.hat-U))
	

